Update: Accurate Determinations of a s from Realistic Lattice QCD 
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We use lattice QCD simulations, with MILC configurations (including vacuum polarization from 
u, d, and s quarks), to update our previous determinations of the QCD coupling constant. Our 
new analysis uses results from 6 different lattice spacings and 12 different combinations of sea-quark 
masses to significantly reduce our previous errors. We also correct for finite-lattice-spacing errors 
in the scale setting, and for nonperturbative chiral corrections to the 22 short-distance quantities 
from which we extract the coupling. Our final result is ay {7 -5 GeV,n/ = 3) = 0.2120 (28), which is 
equivalent to a^g(Mz, Tlf =5) = 0.1183 (8). We compare this with our previous result from Wilson 
loops, which differs by one standard deviation. 
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I. INTRODUCTION 

An accurate value for the coupling constant a s in quan- 
tum chromodynamics (QCD) is important both for QCD 
phenomenology, and as an input for possible theories be- 
yond the Standard Model. Some of the most accurate val- 
ues for the coupling constant come from numerical simu- 
lations of QCD using lattice techniques, when combined 
with very accurate experimental data for hadron masses. 
In this paper we update our previous determinations of 
the coupling from Wilson loops in lattice QCD [l]. Our 
new analysis takes advantage of new simulation results, 
from the MILC collaboration, that employ smaller lattice 
spacings a. We also now account systematically for chiral 
corrections associated with the masses of sea quarks in 
the simulation, and for 0(a n ) uncertainties in the values 
we use for the lattice spacing. 

Few-percent accurate QCD simulations have only be- 
come possible in the last few years, with the development 
of much more efficient techniques for simulating the sea 
quarks; see, for example, [2| for an overview and refer- 
ences. The simulations we use include only light quarks 
(u, d and s) in the vacuum polarization; the effects of c 
and b quarks are incorporated using perturbation theory, 
which is possible because of their large masses. Our lat- 
tice QCD analysis proceeds in two steps. First the QCD 
parameters — the bare coupling constant and bare quark 
masses in the Lagrangian — must be tuned. For each 
value of the bare coupling, we set the lattice spacing to 
reproduce the correct T'-T meson mass difference in the 
simulations, while we tune the u/d, s, c and b masses 
to give correct values for m 2 , 2m 2 K — m 2 , m Vc , and mx , 
respectively; more information can be found in Q. For 



efficiency we set m u — to^; this leads to negligible errors 
in the analysis presented here. Once these parameters 
are set, there are no further physics parameters, and the 
simulation will accurately reproduce QCD. 

Having an accurately tuned simulation of QCD, we 
use it to compute nonperturbative values for a variety of 
short-distance quantities, each of which has a perturba- 
tive expansion of the form 



Y = J2- 



n«y 



{d/a) 



(1) 
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where c„ and d are dimensionless a-independent con- 
stants, and ay(d/a) is the (running) QCD coupling con- 
stant, with rif = 3 light-quark flavors, in the V scheme [3, 
4]. Given the coefficients c n , which are computed using 
Feynman diagrams, we choose ay id /a) so that the per- 
turbative formula for Y reproduces the nonperturbative 
value given by the simulation. Given d and a, and the 
c and b masses, we can then use perturbation theory to 
convert avid/ a) to the more conventional coupling con- 
stant a-^a(Mz,nf = 5), evaluated at the mass of the 
Z meson [5|, Q . 

This analysis is complicated by nonperturbative contri- 
butions to Y and by simulation uncertainties in the value 
of the lattice spacing a, which enters Eq. (fT|). It is also 
complicated by perturbative uncertainties. We know the 
values of the coefficients c„ through order n = 3 (next-to- 
next-to-leading order) for the quantities we examine, yet 
unknown higher-order coefficients still have an impact at 
the level of accuracy we seek. A main focus of this paper 
is to address these complications, and quantify the uncer- 
tainties in our determination of the coupling constant. In 
Section [IT] we review the perturbative expansions for our 
short-distance quantities,all but one of which are derived 
from small Wilson loops [3] . The Monte Carlo simulation 
results for these loops are presented in Section [HTJ We 



discuss finite-lattice-spacing errors and chiral corrections 
in Section ITVl In Section [V] we describe how we combine 
perturbation theory with simulation results using con- 
strained (Bayesian) fitting methods. There we present 
our results and discuss in detail the various uncertain- 
ties that arise. Finally, in Section fVH we summarize our 
results. 



II. PERTURBATION THEORY 

The simplest short-distance quantities to simulate are 
vacuum expectation values of Wilson loop operators: 



W = - 

v ' ran — q 



ReTrPe 



— ig § A-dx 



10), 



(2) 



where P denotes path ordering, A^ is the QCD vector 
potential, and the integral is over a closed ma x na rect- 
angular path. Wilson loops should be calculable in (lat- 
tice QCD) perturbation theory when ma and na are 
small. We computed perturbative coefficients through 
order n — 3 for six small, rectangular loops, and also for 
two non-planar paths: 



BR 



CC 



(3) 



The coefficients for our various loops are derived in [8|. 
The results are for the gluon and quark actions used to 
create the MILC gluon-configuration sets used in this 
study. They also assume n/ = 3 massless sea quarks. 
The quarks in our simulations are not exactly massless, 
but the masses are sufficiently small that the difference 
is negligible, 0(a v (am) 2 ), in perturbation theory (but 
less so nonperturbatively, as we will discuss). 

Perturbation theory is more convergent for the loga- 
rithm of a Wilson loop than it is for the loop itself. This 
is because the perturbative expansion of a loop is dom- 
inated by a self-energy contribution that is proportional 
to the length of the loop, and this contribution expo- 
nentiates for large loops. The length of the loop factors 
out of the expansion when we take the logarithm. This 
structure is evident in Table U where we tabulate the 
perturbative coefficients for the logarithms of our loops. 
The renormalization scales d/a for each quantity are de- 
termined using the procedures described in [3j, \% [9( . 

The perturbative coefficients in \og(W), while greatly 
reduced by the logarithm, are still rather large. They 
can be further reduced in two ways. One is to "tadpole 
improve" W mn by dividing by u where [3| 



u 



(W n ) 1/4 . 



(4) 



The other is to examine Creutz ratios of the loops rather 
than the loops themselves Q. Each procedure signifi- 
cantly reduces the known high-order coefficients, as is 
clear in Table HI We use seven tadpole-improved loops 
and six Creutz ratios in our analysis. Each has smaller 
a v coefficients, which improves convergence, but each 



also has a significantly smaller scale d/a, which slows 
convergence (since ay (d/a) is larger). 

We also include in Table Q] the perturbative expan- 
sion for the tadpole-improved bare coupling constant, 
cciat/Wii, where aut is the coupling constant that ap- 
pears in the gluon action for a given lattice spacing [3]. 
This is another, independent, short-distance quantity 
from which ay can be determined. 

We used Feynman diagrams to compute perturbative 
coefficients c„ for n < 3. Higher-order coefficients can be 
estimated by simultaneously fitting results from differ- 
ent lattice spacings to the same perturbative formula |l[ . 
This is possible because the coupling ay (d/a) changes 
value with different lattice spacings a: 



day(q) 
dq 2 



-f3 Q a v 



PlCXy 



fcay 



fcay 



(5) 



where the Pi are constants [(J. In this paper, we follow 
our earlier analysis by parameterizing the running cou- 
pling by its value at 7.5 GeV, 



a 



a v (7.5GeV,n f = 3). 



(0) 



Given ag , the coupling at any other scale can be obtained 
by integrating Eq. ([5]) (which we do numerically). 

For the purposes of this paper, we define ay in fourth 
order and beyond so that the evolution equation, Eq. ([2|), 
is exact, with no higher-order terms beyond (3s. This 
definition gives precise meaning to the perturbative co- 
efficients c n for n > 4 that we determine by fitting the 
a-dependence of our short- distance quantities [l0( . 

Our main result is a value for a®. To facilitate com- 
parisons with other analyses, we convert this result to 
the MS scheme J6J, add in c and b vacuum polarization 
perturbatively [5j, and then evolve to the mass of the 
Z meson, again using perturbation theory Q. 



III. QCD SIMULATIONS 

The gluon-configuration sets we use were created by 
the MILC collaboration [TlJ. The relevant simulation 
parameters are listed in Table HT1 

The input parameters for a QCD simulation are the 
bare coupling constant and bare quark masses. The cou- 
pling constant is specified through the (3 parameter, listed 
in Table [LTI where 



aiat 



5 
2tt/3' 



(7) 



The bare quark masses, moi(a) for u/d quarks and 
mo s (a) for s quarks, used in the simulations are also 
listed, in units of the lattice spacing and, following MILC 
conventions, multiplied by uq (Eq. ^). The bare masses 
corresponding to fixed physical masses (of, for example, 
pions) vary with the lattice spacing. To facilitate com- 
parisons between lattice spacings, we use first-order per- 
turbation theory to evolve all of our masses to a common 



TABLE I: Perturbative scale and coefficients for several small Wilson loops Wij, Creutz ratios, tadpole-improved Wilson loops, 
and the tadpole-improved bare coupling asiat/Wn. Parameters d and Ci are defined in Eq. |T|. Coefficients Cx,C2,Cz are from 
lattice perturbation theory; coefficients C4, C5 are from the fits to results from multiple lattice spacings described in this paper. 
These results are for the a 2 -improved gluon action used by the MILC collaboration, with the ASQTAD action for vacuum 
polarization from n/ = 3 massless quarks. Similar types of short-distance quantity are grouped. 





d 


Cl 


C2/C1 


C3/C1 


C4/C1 C5/C1 


- log W11 


3.325 


3.06840 


-1.0683(2) 


1.70 (4) 


-4(2) 


-0(4) 


- log W12 


2.998 


5.55120 


-0.8585 (4) 


1.72(4) 


-4(2) 


-1(4) 


- log Wbk 


3.22f 


4.83425 


-0.8547(3) 


1.80 (4) 


-4(2) 


-1(4) 


— log Wcc 


3.047 


5.29758 


-0.7941 (3) 


1.86 (4) 


-4(2) 


-1(5) 


- log Wn 


2.934 


7.87656 


-0.7437(8) 


1.75(5) 


-4(2) 


-1(4) 


- log W14 


2.895 


10.17158 


-0.6870 (8) 


1.70(6) 


-4(2) 


-1(4) 


— log W22 


2.582 


9.19970 


-0.6923(10) 


1.86(5) 


-4(2) 


-1(4) 


— log W23 


2.481 


12.34282 


-0.5995 (13) 


2.00(6) 


-4(2) 


-1(5) 


-logW / l 3 /W 2 2 


2.397 


-1.32313 


0.5969 (84) 


1.11(21) 


-2(2) 


-1(3) 


-\ogW^W 2 2/Wi 2 


2.169 


1.16569 


0.7361 (86) 


1.21(22) 


-2(2) 


-1(3) 


-log Wcc WW W^ 


2.728 


0.92665 


2.2825 (19) 


0.78 (9) 


-4(4) 


-2(6) 


-log Wcc/Wbb.. 


2.730 


0.46333 


0.5103 (35) 


1.16(12) 


-2(2) 


-1(3) 


-\ogWu/W23 


2.066 


-2.17124 


0.5838 (84) 


1.83(29) 


-3(3) 


-1(4) 


-\0%WllW23/Wl2Wvi 


1.970 


1.98345 


0.7062 (88) 


1.64(27) 


-3(3) 


-1(4) 


-logW^M 


2.470 


0.94861 


0.6011(19) 


0.05(8) 


-3(2) 


-1(2) 


-logWBR/'Uo 


2.720 


0.23166 


4.0516(41) 


0.36(16) 


-8(6) 


-3(10) 


-log Wcc /uo 


2.730 


0.69499 


1.6925 (20) 


0.91(8) 


-3(3) 


-1(4) 


-logWia/uo 


1.888 


1.73977 


0.4019 (34) 


-0.44(10) 


-2(1) 


-1(2) 


- log Wu/ul 


1.892 


2.50059 


0.4817(33) 


-0.68(15) 


-2(1) 


-1(2) 


-\ogW 2 2/Uo 


2.290 


3.06291 


0.6149(30) 


0.44 (9) 


-2(2) 


-1(2) 


-logW 23 /«o 


2.030 


4.67183 


0.5714(35) 


0.55(11) 


-2(2) 


-1(2) 



aiat/Wi: 



3.325 



1.00000 



-0.4212(2) 



0.72 (4) 



-4(1) 



-1(2) 



value for the lattice spacing, which we take to be the 
smallest lattice spacing in our analysis: 



m q = mo q (a m [n) 



(8) 



The s-quark masses here are approximately correct. The 
u/d masses are generally too large, but small enough to 
allow accurate extrapolations to the correct values. 

The lattice spacing is not an input to QCD simula- 
tions. Rather it is extracted from calculations of physical 
quantities in the simulation. Here we use MILC's deter- 
minations of r\/a for this purpose, where r\ is defined in 
terms of the static-quark potential [ll(. The values for 
each configuration set are listed in Table [Til To obtain 
the lattice spacing, we need to know r\. We use the value, 
r\ = 0.321 (5)fm, determined from simulation results for 
the T'-T mass splitting [1^ . The uncertainties quoted 
for r\ja in Table HTl are predominantly statistical; they 
do not include potential errors due to the finite lattice 
spacing or mistuned light-quark masses, which we will 
discuss later. 

The lattices we use here have lattice spacings that 
range from 0.18 fm to 0.045 fm. The spatial volumes are 
2.4 fm across or larger in each case. 

Our simulation results for the vacuum expectations of 
our 8 different Wilson loops, each for each of our 12 dif- 
ferent configuration sets, are presented in Table ITTT1 The 



TABLE If: QCD parameters for the 12 different sets of gluon 
configurations used in this paper [Tj. Parameter /3 specifies 
the bare coupling constant. The inverse lattice spacing is 
specified in terms of the r\ , and the bare quark masses are in 
units of the lattice spacing, and multiplied by Uq. The spatial 
and temporal sizes, L and T, axe also given. Configuration 
sets that were tuned to have the same lattice spacing are 
grouped. 



Set 


P 


n/a 


auoinoi 


auomos 


L/a 


T/a 


1 


6.458 


1.802(10) 


0.0082 


0.082 


16 


48 


2 


6.572 


2.133(14) 


0.0097 


0.0484 


16 


48 


3 


6.586 


2.129(12) 


0.0194 


0.0484 


16 


48 


4 


6.76 


2.632(13) 


0.005 


0.05 


24 


64 


5 


6.76 


2.610(12) 


0.01 


0.05 


20 


64 


6 


6.79 


2.650(08) 


0.02 


0.05 


20 


64 


7 


7.09 


3.684(12) 


0.0062 


0.031 


28 


96 


8 


7.11 


3.711(13) 


0.0124 


0.031 


28 


96 


9 


7.46 


5.264(13) 


0.0018 


0.018 


64 


144 


10 


7.47 


5.277(16) 


0.0036 


0.018 


48 


144 


11 


7.48 


5.262(22) 


0.0072 


0.018 


48 


144 


12 


7.81 


7.127(34) 


0.0028 


0.014 


64 


192 



uncertainties quoted are statistical. Step-size errors, due 
to the algorithm used to generate gluon configurations, 
are no larger than the statistical errors [l3| and therefore, 
like statistical errors, are negligible; we will ignore them 
here. 



IV. SYSTEMATIC ERRORS 

The goal of our analysis is to determine u$ = 
av(7-5GeV). The only relevant systematic errors, other 
than from the truncation of perturbation theory, are from 
nonperturbative effects and from a 2 errors in our deter- 
mination of the lattice spacings. Finite-volume errors 
are no larger than our statistical errors, as we have ver- 
ified by examining configuration set 5 with L/a = 28 in 
addition to L/a — 20. Statistical errors are also neg- 
ligible (and therefore we ignored statistical correlations 
between different Wilson loops when computing Creutz 
ratios, whose real statistical errors are 2-3 times smaller 
than what we use here). We consider each systematic 
effect in term. 



A. Chiral Corrections 

Wilson loops, being very short-distance, are almost in- 
dependent of the light-quark masses. The dependence 
in perturbation theory is 0{ay(am q ) 2 ), which is negli- 
gible here given other errors. There is a larger contribu- 
tion, however, from nonperturbative contributions that 
is important to our analysis. This contribution can be 
parameterized using chiral perturbation theory and the 
operator product expansion, which says that an arbitrary 
QCD operator Oqcd that is local at scale A can be ex- 
panded in terms of local operators O n from the chiral 
theory: 



and U = exp(i(f>/F) with 



O 



QCD 






On 

■■ A d n 



(9) 



where d n is the dimension of O n minus the dimension of 
Oqcd • Here equivalence between the left-hand and right- 
hand sides means that matrix elements of the operators 
are equal for comparable physical states in QCD and the 
chiral theory. 

For Wilson loops, we are interested in vacuum expec- 
tation values and singlet operators. The scale A for a 
loop of size L is A ~ 1/L. Consequently we expect 



W = b + b x L Tr (m(U + [/ f )) 
+ 6 2 £ 2 Tr(<9 A1 l^[/ t )- 



(10) 



r T °/V2- 



Vs/VQ 



K' 



7T T 
-7T%/2 4 

K° 



WV6 



K+ ' 
K° 

-2VV6. 

(11) 



and F w 92 MeV. 

Taking the vacuum expectation value and a logarithm, 
and keeping only the leading 0(a) terms, we get 

log(W) » w (0) (l + w£>a (Tr (m(U + [/+)))' 

« w {0) (l + w\l ) a{2mi + m s ) + ■ ■ ■ ) . (12) 

Standard methods can be used to compute higher-order 
corrections, including chiral logarithms, from the expan- 
sion of Tr (m(U + U^)), but these are too small to be 
relevant to our analysis. 

The leading contribution, w^°\ is obtained from the 
perturbative analysis discussed in Section HH provided 
the loops are sufficiently small to be perturbative. We 
expect Wm to be roughly independent of loop size since 
«/ ' is approximately proportional to L/a (see Sec- 
tion HT}. 

We can estimate the size of Wm from a simple argu- 
ment. For light-quark hadrons, hadronic quantities like 
meson decay constants or baryon masses depend approx- 
imately linearly on the masses of their valence quarks. 
The mass m v of a valence quark makes a contribution of 
order Q m v /A to some hadronic quantity Q, where A is a 
momentum scale characteristic of the size of the hadron 
(« the chiral scale, for light-quark hadrons). From ratios 
of decay constants like Jk/ fiv or of baryon masses like 
m(A°)/m(p + ), it is clear that m s /A is of order 20%, and 
therefore that A w 400 MeV. Empirically contributions 
from individual sea-quark masses are 3-5 times smaller 
than those from individual valence-quark masses [14j . 
Consequently the relative contribution from a sea-quark 
mass m q should be roughly m.j/1.2 GeV. 

Now consider Wilson loops. The m q dependence of 
log(Wn), for example, should be much smaller than that 
for a light-quark hadron because the loop is much smaller 
than the hadron. The typical radius of such hadrons is 
around 1 fm, so we expect the relative contribution to 
Wn from a sea-quark mass of m q to be approximately 



lfm 1.2 GeV 



am q 
~6~ 



(13) 



where m = dia,g(m u ,md,m s ) breaks chiral symmetry, 



Therefore we expect Wm' — 0(1/6). This implies correc- 
tions to our log(W)s, for example, of order 1-2% on the 
coarsest lattices and 0.3-0.5% on the finest lattices — 
which is large compared with the statistical errors in 
these quantities, and therefore important. 

In most lattice calculations we want the light-quark 
masses as close to their physical values as possible, so 



TABLE III: Simulation results for the vacuum expectation values of various small Wilson loops. Results are given for each of 
the 12 different configuration sets in Table HT1 



Set Wu 



Wv 



Wv. 



Wv 



W-, 



W 2 ; 



w B 



Wcc 



1 


0.534101(17) 


0.280720(22) 


0.149263(21) 


0.079710(19) 


0.087438(23) 


0.030150(16) 


0.338982(22) 


0.287376(25) 


2 


0.548012(51) 


0.298624(68) 


0.165063(67) 


0.091701(63) 


0.101572(73) 


0.038333(54) 


0.356763(68) 


0.306315(78) 


3 


0.549470(53) 


0.300310(70) 


0.166530(70) 


0.092797(63) 


0.102640(76) 


0.039007(54) 


0.358570(70) 


0.308140(79) 


4 


0.567069(16) 


0.323163(21) 


0.187281(24) 


0.109122(27) 


0.121542(19) 


0.050751(15) 


0.381148(25) 


0.332184(27) 


5 


0.566961(21) 


0.322987(27) 


0.187084(26) 


0.108927(21) 


0.121341(29) 


0.050579(21) 


0.380988(26) 


0.331996(29) 


6 


0.569716(21) 


0.326496(27) 


0.190278(26) 


0.111479(21) 


0.124204(29) 


0.052397(21) 


0.384479(26) 


0.335685(29) 


7 


0.594843(7) 


0.359761(9) 


0.221624(10) 


0.137271(10) 


0.153433(12) 


0.072261(10) 


0.417002(9) 


0.370239(10) 


8 


0.596408(12) 


0.361838(19) 


0.223616(17) 


0.138946(16) 


0.155315(18) 


0.073593(16) 


0.419020(16) 


0.372372(17) 


9 


0.620813(5) 


0.394837(8) 


0.255897(9) 


0.166723(9) 


0.186116(7) 


0.096208(6) 


0.450947(7) 


0.406300(8) 


10 


0.621462(3) 


0.395717(4) 


0.256770(5) 


0.167486(5) 


0.186959(5) 


0.096852(5) 


0.451798(4) 


0.407210(5) 


11 


0.622115(2) 


0.396607(4) 


0.257650(4) 


0.168257(4) 


0.187809(5) 


0.097491(4) 


0.452654(3) 


0.408123(4) 


12 


0.641947(2) 


0.423992(3) 


0.285304(5) 


0.192943(5) 


0.214759(4) 


0.118532(3) 


0.478903(3) 


0.436064(4) 



that lattice results reproduce what is seen in experi- 
ments. The situation for our Wilson loops is different, 
however. In our simulations here we are trying to isolate 
the perturbative part of the loop, in order to compare 
it with perturbation theory (not experiment), and the 
linear quark-mass dependence is a nonperturbative con- 
tamination that we want to remove. Consequently the 
precise values of the quark masses are not relevant so 
long as they are small enough that we can correct for 
them (or ignore them), which is the case here. 



While a leading-order condensate value of 0.006, for ex- 
ample, shifts logW23 by about 25% for our largest lat- 
tice spacings, the shift is less than 0.1% for the smallest 
lattice spacing, which is more important in our analy- 
sis. Smaller loops are much less sensitive: for example, 
this gluon condensate shifts log Wu by only 0.3% for the 
largest lattice spacings, and by only 0.003% for the small- 
est lattice spacings. The two Creutz ratios that involve 
W23 are the most sensitive to condensate contributions, 
but even they are shifted by only 0.2-0.25% for the small- 
est lattice spacings J17I ]. 



Gluon Condensate 



The leading gluonic nonperturbative contribution 
comes from the gluonic condensate, (a s G 2 /tt). The con- 
tribution of the condensate to a Wilson loop is easily 
calculated to leading order in perturbation theory: 

SW cond = -^ (£\ a\a s G 2 /v) (14) 

where A is the loop area for planar loops. We remove 
this contribution from our Wilson loops before compar- 
ing them with perturbation theory. The value of the 
condensate is not well known, so we take (a s G 2 /it) = 
0.0 ± 0.012 GeV 4 , which covers the range of expecta- 
tions [lH . We also allow for higher-dimension condensate 
contributions by replacing 



SW a 



SW cwd {l + w£ nd (aA g ) 



™z d (^ a ) 4 



(15) 



where we take A c 



.,(0 



1 GeV and coefficients «^ o ; nd = Oil. 



To be certain that we do not underestimate errors we 
include 10 condensate terms in all [16]. 

We chose the number of condensate terms here some- 
what arbitrarily. Only results from the largest loops are 
affected appreciably even by the leading-order conden- 
sate correction, and then only by amounts of order a 
standard deviation in our final results for the coupling. 



C. Finite-a Errors 

In our analysis, the scale for the couplings comes from 
the lattice spacing, and the lattice spacing comes from 
measurements of r\/a in the simulations. As for any 
physical quantity, lattice QCD measurements of n have 
finite-a errors; and, using an analysis similar to the one 
we outlined for Wilson loops, they should also be approx- 
imately linear in the sea-quark masses. Consequently we 
expect 



„lat 



where: 



„( 2 ) 



„(i) 



r 1 (l + r^(a/r i y+r\^r 1 (2Sm l + Sm s ) + - ■ • ) (16) 



(2) 
1 

la 



0(a s fa 1/3) [18| . since the gluon action 
has no tree-level errors in 0(a 2 ); and r lr ^ = 0(1/6), 
following the discussion for Wilson loops. Here 5m q is the 
simulation's tuning error in the mass for sea-quark q — 
dmi ss mi for our simulations, while Sm s sw 0. These 
corrections could affect our lattice spacings by as much as 
several percent, although the impact on ao is suppressed 
by a power of ao and so is much less. We allow for both 
corrections in our analysis. 



V. ANALYSIS AND RESULTS 

We have 22 different short-distance quantities in our 
analysis, each of which produces a separate value for 



a = ay (7.5 GeV). These consist of \og(W)s for each of 
8 Wilson loops, 6 independent Creutz ratios built from 
these loops, 7 tadpole-improved \og(W)s, and the tadpole 
improved bare coupling ai a t/Wn. We have 12 values for 
each of these quantities, with one for each configuration 
set in Table HU In this section we discuss first the fitting 
method used for extracting ao, and then we review our 
results. 



A. Constrained Fits 

We analyze each short-distance quantity Y separately. 
We use a constrained fitting procedure, based upon 
Bayesian ideas 19], to fit the values Yj ± <jy i coming 
from each of our configuration sets (Table IIIIj) to a single 
formula. In this procedure we minimize an augmented 
X 2 function of the form 

2 _ y^ (Y -Y(a l ,(am q )i,a 0l y { r n l ,Cn 1 d)) 2 .\r^ x2 



i=i 



a 



Yi 



where i labels the configuration set, and 
Y(a i ,(am q )i,aQ, y£\c n , d) 



(17) 
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= \l + y { ^{2arni + am s )i\^c n av{d/ai). (18) 



71=1 



The sea-quark mass dependence here is from Eq. (|12|l. 
The lattice spacing in each case is determined from the 
simulation values for {r\/a)i from each configuration set 
(Table HH using 



r\ 



(n/o) 



- (l + ri^o/n)? + r^ (2nm0i) , (19) 



which follows from Eq. p6[) . taking <5m; sa to/ and 
Sm s w 0, and r\ = 0.321(5) fm [12]. Here (rim q )i = 
(am q )i{ri/a)i. Given the lattice spacing, the coupling 
ay (d/a) is computed from ao by integrating Eq. ([5]) nu- 
merically. 

The \ 2 function is minimized by varying fit parameters 
like the c n (but not d which is effectively exact). Every 
fit parameter in our procedure is constrained by an extra 
term or "prior" 8x\ in the x 2 function. The expansion 
parameters c n from perturbation theory, for example, are 
constrained by 



5x1 



10 

E 

ra=l 



(c„ 



which implies that the fit will explore values for c n that 
are centered around c„ with a range specified by a c „ '■ c„± 
cr Cn . For n < 3, we set c n to the value obtained from our 
numerical evaluation of the relevant Feynman diagrams, 
with a Cn equal to the uncertainty in that evaluation. For 
n > 4, we set c„ = and 



2.5max(|ci|,|c 2 |,|c 3 | 



(21) 



Thus the c„s in the fit are constrained by the values ob- 
tained from our Feynman integrals where these are avail- 
able (taking correct account of the uncertainties in those 
values) , while the others are allowed to vary over a range 
that is 2.5 times larger than the largest known coefficient. 
The factor 2.5 was chosen using the empirical Bayes cri- 
terion, described in |19j, applied to the log(W)s; apply- 
ing the same criterion to the other quantities would have 
given smaller factors, but we take the more conservative 
factor of 2.5 for these as well. 

We include seven c n s beyond the ones currently known 
from perturbation theory to illustrate an important issue. 
In reality there are infinitely many c„s, but in practice 
the various uncertainties in our analysis mean that it is 
sensitive only to the first few. As we add c„s the fit im- 
proves but only up to a point — n — 4 for \og{W)s. As 
long as priors are included in x 2 ■, terms can be added 
beyond this point but they have no effect on the result 
of the fit (including the error estimate) or on the quality 
of the fit. We add terms through n = 10 to be certain 
we have reached this point. Our analysis is not suffi- 
ciently accurate to yield new information about c„s with 
n > 4 (beyond what is incorporated in the prior); but, by 
adding enough c„s so that the fit results and errors cease 
changing, we guarantee that our final error estimates in- 
clude the full uncertainty due to the fact that we have a 
priori values for only a few of the coefficients. 

Other fit parameters, like ao, ym , r[ a , and r^, must 
also have priors: 



Sx 2 o 



, _(log(a )-log(a )) 2 



+ 



log(cKo) 

Via 'la) 



(v {1) -v {1) ) 2 
\ym ym ) 



^(^2 



( r v> -r ( >) 

\ lm XmJ 



(22) 



,( 2 ) 



We constrain log(ao) to be —1.6 ± 0.5; this prior has 
negligible effect on the fits because it is so broad (and 
the fits are so sensitive to ao). Following the discussion 
in Section HVl we set 







-(i) = -(i) = 

Urn lm 



v (i) = °v<i> = 1/6- 



(23) 



We checked the width of these two priors using the em- 
pirical Bayes criterion and found that, in fact, this is the 
optimal width indicated by our simulation results. For 
r la _ , the empirical Bayes criterion suggests a width for the 
(20) prior that is twice what we anticipated in Section HV CI 



-(2) 

i 
la 







cr (2) 



2a s w 0.6. 



(24) 



We use this more conservative prior in our fits. Higher- 
order corrections are easily added but have no impact 
because the corrections are too small to matter, given 
the size of our other errors. 

Our simulation result for (ri/a)j, which is used to de- 
termine the lattice spacing a, for the i th configuration set 



(Eq. (|19p ). is not exact. To include its uncertainty in our 
analysis we treat (ri/a), as a fit parameter, to be varied 
while minimizing \ 2 j but with a prior whose mean is the 
value measured in the simulation and whose width is the 
measured uncertainty (as in Table HI)) . We can incorpo- 
rate the uncertainty in the value of r\ using the same 
trick, with T\ as a fit parameter: 



Sxli = 



(ri -ri) 5 



12 

£ 

8 = 1 



((n/ah - (n/a).,) 2 

o 
(ri/a)i 



(25) 



where r x ± cr ri = 0.321 ± 0.005 fm [jj. 

The c and b masses are required to convert ag to 
ajj^(Mz,nf — 5). We account for the uncertainties in 
these masses by including them as fit parameters, with 
appropriate priors, together with fit parameters for un- 
known high-order terms in the MS /3-function, and in the 
perturbative formulas for incorporating c and b vacuum 
polarization |5|, Q. For the /3-function, we allow for a 
sixth-order term /^a^- in the evolution equation (anal- 
ogous to Eq. iJSJl for ay) where /?4 is a fit parameter with 
a prior centered on /3 4 = with width 



*fU=-nBx(\0o\,\Pi\M,\fo\) 



(26) 



for the MS /3jS. We include analogous corrections, fit 
parameters and priors for the formulas for c and b vacuum 
polarization. 



B. Results 

The results from our 22 determinations of the coupling 
are listed and shown in Figure [1] The gray band corre- 
sponds to our final result of 



'MS 



(M z ,n/ = 5) = 0.1183(8) 



(27) 



which was obtained from a weighted average of all of 
22 determinations [2(|. Our error estimate here is that 
of a typical entry in the plot; combining our results does 
not reduce errors because most of the uncertainty in each 
result is systematic. The individual results in the plot are 
consistent with each other: x 2 /22 = 0.2 for the 22 entries 
in Figure [T] And the fits for each quantity separately are 
excellent as well: x 2 /11 = 0.3 to 0.6 for our fits to the 
12 pieces of simulation data (one from each configuration 
set) for each quantity. The results in Figure[T]are derived, 
using perturbation theory (Section[TlJ , from the fit values 
for ao, which average to 

a Q =a y (7.5GeV,n / = 3) = 0.2120(28), (28) 

where again the error is that of a typical result for a single 
short-distance quantity (it is not reduced by one over the 
square root of the number of inputs) . 

Figure [2] reveals more details about our fit. The top 
panel in this figure shows the values of ay (d/a) coming 




_i_ 



_i_ 



0.1185(8) 
0.1185(8) 
0.1185(8) 
0.1184(8) 
0.1184(8) 
0.1184(9) 
0.1183(9) 
0.1181(9) 



logWn 

log W12 

log W B r 
log W C c 
log W 13 

log W u 

log W 22 
log W 23 



0.1186(9) \ogW 13 /W 22 
0.1185(9) logWnW 2 2/Wi2 

0.1185(9) logW CC W BR /W^ 
0.1184(8) \ogW cc /W BR 
0.1171(11) log Wu/W 2 3 
0.1174(11) \ogWuW 23 /W 12 W 13 



0.1183(7) 

0.1187(11) 

0.1184(9) 

0.1181(7) 

0.1186(8) 

0.1183(8) 

0.1176(8) 
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FIG. 1: Values for the 5-flavor Ojjg at the Z-meson mass from 
each of 22 short-distance quantities. The gray band indicates 
our final result, 0.1183 (8). % 2 P er data point is 0.2. 



from every short-distance quantity for every lattice spac- 
ing in our configuration sets. The ays plotted here were 
obtained by refitting each piece of simulation data sep- 
arately, rather than fitting results from all lattice spac- 
ings simultaneously as above. In these fits we used the 
values for c„ with n > 3, win , etc. obtained from our 
simultaneous fit to all lattice spacings [21(, which is why 
the individual data points align well with the perturba- 
tive result for ay (d/a) (the gray band). The fact that 
different points align so well is an indication of the self- 
consistency of our perturbative analysis across all scales 
and for all quantities. The size of the error bars for dif- 
ferent points is determined by the perturbative and non- 
perturbative uncertainties associated with each piece of 
simulation data. Points with error bars much larger than 
the uncertainties in the perturbative ay (that is, much 
larger than the vertical width of the gray band) have lit- 
tle impact on our overall fits. The bulk of the uncertainty 
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FIG. 2: Values for av versus d/a from each short-distance 
quantity at each lattice spacing, with and without corrections 
for gluon condensates. The gray band shows the prediction 
from QCD evolution (Eq. ([5])) assuming our composite fit 
value (Eq. (gSJ). 



at low momentum comes from uncertainties in the gluon 
condensates. This is obvious when the results are rean- 
alyzed without corrections for the condensates (bottom 
panel in Figure [21) ■ The most important simulation data 
is at large d/a, where errors are smaller than the plot 
points whether or not condensates are included. 

It is useful to separate our error estimates into compo- 
nent pieces. The error estimate produced by our fitting 
code for a quantity like cujjjg is approximately linear in 
all the variances a 2 that appear in the \ 2 function: 



12 



10 
n=l 



+ c r m aim + c i2) o 2 j 2) + 



(29) 



This works when errors are small, as they are here. To 
isolate the part of the total error that is associated with 
the statistical uncertainties in the Yi, for example, the fit 
is rerun but with the corresponding variances rescaled by 
a factor / close to one (/ = 1.01, for example): 



for i = 1 ... 12. Then 
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(31) 



4 



f°Yi 



(30) 



The square root of this quantity is the part of the total 
error due to the statistical uncertainties in the Y t . This 
procedure can be repeated for each prior or group of pri- 
ors that contributes to the x 2 function. The sum of the 
variances obtained in this way for each part of the total 

error should equal a 2 a _; if it does not, errors may not 

be sufficiently small to justify the linear approximation 
in Eq. ([29]) [2j. 

In Table [TV] we present error budgets computed in this 
fashion for a sample of our determinations of a-y^{Mz)- 
This table shows that our largest errors come from un- 
certainties in the perturbative coefficients with n > 4, 
statistical errors in the simulation values for (ri/o),, sys- 
tematic uncertainties in the physical value for n, and 
finite-a lattice errors in r±. Uncertainties in the param- 
eters used to convert ao = ay (7.5 GeV, n/ = 3) into 
ajjg(Mz,nf = 5) have negligible impact. Also negligi- 
ble are uncertainties due to the gluon condensate, and 
statistical errors in the Wilson loops. 

Our errors are greatly reduced because we can bound 
the size of perturbative coefficients c n for n — 4 and 
beyond. This is possible because we are fitting simula- 
tion data from six different lattice spacings simultane- 
ously. As noted in PJ, the n = 4 coefficients are large, 
particularly for log(IV)s where typically our fits imply 
c 4 /ci sw —4(2). As expected, perturbative higher-order 
coefficients are smaller for other quantities: for exam- 
ple, we find typically c^/ci ~ —2(2) for tadpole-improved 
loops. The fit results for Gi/ei and C5/C1 for each of our 
short-distance quantities are given in Table [TJ 

We tested the stability of our analysis procedure in 
several ways: 

• Discarding simulation data: Dropping data for 
any one of the lattice spacings gives results that 
are almost identical to our final result: the value 
of a-^g(Mz) varies by no more than 0.12% from 
our final result, and its uncertainty ranges be- 
tween 0.00083 and 0.00093. Dropping the two 
smallest lattice spacings, which are the most im- 
portant, shifts a^g(Mz) to 0.1176(14). Keeping 
just the four, three and two smallest lattice spac- 
ings gives 0.1183(9), 0.1180(10), and 0.1179(10), re- 
spectively (for sets 4-12, 7-12, and 9-12). 

• Perturbation theory scale changes: Our results do 
not depend strongly on the choice of scale d/a used 
in the perturbation theory for each quantity. Re- 
expanding our perturbation theory for d — > d/1.5 or 
d — > 1.5d, for example, shifts the overall a^g(Mz) 
to 0.1181(8) or 0.1184(8), respectively [1|. 

• MS throughout: Re-expressing the perturbation 
theory for each quantity in terms of a^jg in 





log Wii 


log Wl2 


log W22 


logWnW 2 2/W? 2 


logWVo 


lOgW / 22/«0 


Olat/Wll 


C1...C3 


0.1% 


0.1% 


0.1% 


0.3% 


0.1% 


0.1% 


0.1% 


c n for n > 4 


0.2 


0.3 


0.3 


0.4 


0.3 


0.4 


0.3 


aro„ rim g extrapolation 


0.1 


0.1 


0.0 


0.1 


0.1 


0.1 


0.0 


(o/ri) 2 extrapolation 


0.2 


0.3 


0.4 


0.3 


0.2 


0.2 


0.0 


ij'i/a)i errors 


0.4 


0.4 


0.4 


0.3 


0.3 


0.3 


0.3 


ri errors 


0.3 


0.3 


0.3 


0.3 


0.3 


0.3 


0.3 


gluon condensate 


0.1 


0.1 


0.1 


0.2 


0.1 


0.1 


0.1 


statistical errors 


0.0 


0.0 


0.0 


0.1 


0.0 


0.0 


0.0 


V -> MS -*■ M z 


0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


Total 


0.6% 


0.6% 


0.7% 


0.7% 


0.6% 


0.6% 


0.5% 



TABLE IV: Sources of uncertainties in determinations of Ojjb^Mz, n/ = 5) from various short-distance quantities. Uncertainties 
are given as percentages of the final result in each case. 



place of ay gives almost the same overall results, 
0.1185(10), but leads to significantly larger high- 
order coefficients in perturbation theory (2.5 times 
larger for small loops), somewhat greater dispersion 
between results from different quantities (x 2 /22 
of 0.5 instead of 0.2), and larger uncertainties in 
the results from most quantities. The scale-setting 
procedure used to select the ds is tailored specifi- 
cally for ay expansions; this is reflected by these 
results. 



• Adding more/fewer perturbative terms: We allow 
terms up through tenth order in the perturbative 
expansions for the various short-distance quanti- 
ties. Adding further terms has no impact on our 
results. Restricting perturbation theory to only 
fourth or fifth order also leaves our final result un- 
changed. Fitting is impossible with fewer than four 
terms: with three terms fits for individual Wilson 
loops, for example, to data from all 12 configura- 
tion sets are poor, with x 2 /12 becoming as large 
as 1.9 (rather than 0.4); and the couplings coming 
from the 22 different short-distance quantities dis- 
agree with each other, giving x 2 /22 = 1.45 (rather 
than 0.16). 



Adding more/fewer nonperturbative terms: Adding 
higher-order terms in the chiral expansion in sea- 
quark masses (Eq. 1]12[) ) or further terms in the 
gluon-condensate expansion (Eq. ([15])) does not 
change our final result at all. Omitting all correc- 
tions for the gluon condensates increases a.y^{Mz) 
by two thirds of a standard deviation, to 0.1189(7). 
If we keep only the three smallest lattice spacings, 
which are the least sensitive to nonperturbative ef- 
fects, we get 0.1180(10) whether or not the gluon 
condensates are included. We cannot fit all of our 
simulation data if we omit the chiral correction. 
Fitting without chiral corrections becomes possi- 
ble if we keep only the subset of our data with 
m u/d/ m s ~ 0.2 (sets 2, 5, 7, 10, and 12); this gives 



a m (M z ) = 0.1181(9). (Our fit to log(Wn) gives 



u# = -0.18(6) r £> =-0.08(8), (32) 

which is typical of the other fits.) 

Each of the variations examined here gives results that 
agree with our final result to within a standard deviation, 
suggesting that we have not underestimated the uncer- 
tainty in our result. 

Our new result is one standard deviation above our 
previous result from Wilson loops [l(, a^g(Mz) = 
0.1170(12), and has an error that is 33% smaller. Our 
new analysis differs in two important ways from our ear- 
lier work. First we include more lattice spacings, includ- 
ing one that is 50% smaller than the smallest we used 
before. (We used only configuration sets 1, 5 and 7 be- 
fore.) This significantly reduces the errors. Second we 
now use more accurate values for r\/a. These reduce un- 
certainties in the ratios of lattice spacings from different 
configuration sets, to a third of what they were in our 
earlier analysis. This matters since comparing results at 
different lattice spacings bounds the uncalculated high- 
order perturbation theory coefficients in our analysis (c n 
for n > 4) . We are also allowing for larger finite-a errors 
in r±/a on the coarsest lattices than we did previously. 
The changes in ri/a, together with the smaller lattice 
spacing, account for most of the increase in our final re- 
sult. 

Another change, which has less impact, is the inclusion 
of possible higher-dimension condensates. We also now 
do a more systematic analysis of effects due to the sea- 
quark mass, fitting results with many different masses, 
but the effect on our final result is small. Finally, we now 
use better scales d/a for the Creutz ratios and tadpole- 
improved loops than in our previous analysis [9j. Using 
the new scales shifts our final result up by only a third of 
a standard deviation, but the dispersion between results 
from different short-distance quantities is decreased from 
X 2 /22 = 0.6 to 0.2. 
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VI. CONCLUSIONS 

Any high- precision determination of a s based upon lat- 
tice QCD simulations has to address several key issues: 

• Finite-Lattice-Spacing Errors: Errors due to the fi- 
nite lattice-spacing can enter in two ways. First 
they affect lattice determinations of the physical 
quantity or quantities used to set the scale of the 
coupling. In our analysis we use simulation val- 
ues for ri/a, from the static-quark potential, to 
determine ratios of scales from different configu- 
ration sets, and simulation values for the T'-T 
mass difference to set the overall scale (12J , In 
each case we use data from multiple lattice spac- 
ings to bound finite-a errors, which are small be- 
cause we use highly-improved discretizations in 
our simulations. The second source of finite-a er- 
rors, for some analyses (but not ours), is the lat- 
tice determination of the short-distance quantity 
that is compared with perturbation theory (to ex- 
tract a s ). A short-distance quantity that is de- 
fined in continuum QCD — for example, changes 
V'(tv) — V(r;,) in the static-quark potential for small 
rs [ll lS| , or current-current correlators for c-quark 
currents [24j — will have finite-a errors that must 
be included in the final error analysis. The use of 
multiple lattice spacings is again important. This 
is not an issue for us here because we analyze our 
short-distance quantities using lattice QCD pertur- 
bation theory, which treats finite-a effects exactly 
(that is, to all orders in a, order-by-order in ay)- 
Both the simulation results and the perturbation 
theory for our 22 short-distance quantities are free 
of finite-a errors. This greatly facilitates our use 
of results from multiple lattice spacings to bound 
uncalculated higher-order terms from perturbation 
theory. 

• Truncation Errors from Perturbation Theory: The 
coupling is determined by comparing perturbation 
theory with (nonperturbative) simulation results 
for a short-distance quantity. Generally the per- 
turbation theory is known through only a few low 
orders in a s . The error analysis for any determi- 
nation of the coupling must account for the uncal- 
culated (but certainly present) terms from higher- 
order perturbation theory. We not only account 
for the possibility of higher-order terms (through 
tenth order), using our Bayesian priors, but also 
attempt to estimate the size of these corrections 
by comparing values of our short- distance quan- 
tities at five different momentum scales d/a, cor- 
responding to our five lattice spacings. We find 
sizable contributions from high-order terms, par- 
ticularly for \og(W)s: leaving them out would shift 
our final result for the coupling down by one to two 
standard deviations (and lead to poor fits for most 
of our short-distance quantities). The agreement 



between our 22 different short-distance quantities, 
some with very different perturbative expansions 
(see Section|n|, is important evidence that we have 
analyzed truncation errors correctly. 

• Sea-Quark Vacuum Polarization: In our previous 
analysis [l(, we showed that the coupling is quite 
sensitive to contributions from the vacuum polar- 
ization of sea quarks: ay^(Mz) is 30% smaller 
when all quark vacuum polarization is omitted. It 
is therefore important to include vacuum polariza- 
tion from all three light quarks. Vacuum polar- 
ization corrections from heavy quarks (c, 6 and t) 
can be computed using perturbation theory, but 
light quarks (m, d and s) can only be incorporated 
nonperturbatively. In the past we have used simu- 
lations with fewer than three light-quarks and ex- 
trapolated to nf — 3 (l/a^jg(Mz) appears to be 
reasonably linear in n/) [251 ]. Here (and in our ear- 
lier paper Jl|) contributions from all three light- 
quarks are included in the configurations provided 
to us by the MILC collaboration. We also account 
for the small but (barely) measurable dependence 
upon the sea-quark masses. 

• Other Lattice and Nonperturbative Artifacts: Usu- 
ally one must worry about the finite volume of the 
lattice in a QCD simulation. Our Wilson loops, 
however, are about as ultraviolet singular as is pos- 
sible on a lattice, and so are completely insensi- 
tive to the volumes of our lattices (2.5 fm across). 
Another issue, for continuum as well as lattice de- 
terminations of the coupling, is the possibility of 
nonperturbative contributions to the short-distance 
quantity. Our quantities are sufficiently short- 
distance that we do not expect appreciable nonper- 
turbative contamination. We nevertheless allowed 
for nonperturbative contributions from both gluons 
and quarks. The expected size of nonperturbative 
contributions varies widely over our set of 22 differ- 
ent short-distance quantities and 6 different lattice 
spacings. The excellent agreement among all of our 
results is strong evidence that we understand these 
systematic errors. 

In this (and our previous) paper, we have addressed 
all of these issues. We have extended our earlier anal- 
ysis of the strong coupling constant from Wilson loops 
in lattice QCD (and hadronic spectroscopy) to include 
results from 22 different short-distance quantities com- 
puted on 12 different lattices, with 6 distinct lattice 
spacings and a variety of sea-quark masses. We ex- 
tracted a new value for the QCD coupling by compar- 
ing these 22 x 12 = 264 different pieces of simulation 
data, varying by a factor of seven in momentum scales 
(d/a from 2.1 to 14.7 GeV), with perturbation theory. 
Our result, ay^(M z , nf — 5) — 0.1183(8), is in excellent 
agreement with our previous result from Wilson loops |l[ , 
0.1170 (12), and also with non-lattice determinations: for 



11 



example, the world averages 0.1176(20) from [26| and 
0.1189(10) from [23]. Our new result also agrees well 
with our very recent result, 0.1174(12), from current- 
current correlators computed using lattice QCD 24]. 

While they are derived from the Wilson loops, 
our Creutz ratios and tadpole-improved loops provide 
coupling-constant information that is independent from 
that coming from the loops directly. This is because the 
highly ultraviolet contributions that dominate the loops 
largely cancel in the other quantities, making the latter 
more infrared. Consequently both perturbative and non- 
perturbative behavior differs significantly from quantity 
to quantity. This is particularly true of the sensitivity to 
nonperturbative contributions: for example, our most in- 
frared Creutz ratios are more than 100 times more sensi- 
tive to gluon condensates than our most ultraviolet loops. 
That all of our quantities agree on the coupling (Figure[T]) 
is strong evidence that we understand the systematic er- 
rors involved. 

The close agreement of our results with non-lattice de- 
terminations of the coupling is a compelling quantitative 
demonstration that the perturbative QCD of jets, and 
the QCD of lattice simulations, which encompass both 



perturbative and nonperturbative phenomena, are the 
same theory. It is also further evidence that the sim- 
ulation methods we use are valid. While early concerns 
about the light-quark discretization used here have been 
largely addressed [28|, |29[, it remains important to test 
the simulation technology of lattice QCD at increasing 
levels of precision given the critical importance of lattice 
results for phenomenology [301 ]. 
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